Momentum transferred to a trapped Bose-Einstein condensate 
by stimulated light scattering 



A. Brunello 1 , F. Dalfovo 2 , L. Pitaevskii 1,3 , S. Stringari 1 and F. Zambelli 1 

1 Dipartimento di Fisica, Universita di Trento, and 
Istituto Nazionale per la Fisica della Materia, 1-38050 Povo, Italy 
2 Dipartimento di Matematica e Fisica, Universita Cattolica del Sacro Cuore, and 
Istituto Nazionale per la Fisica della Materia, gruppo collegato di Brescia, 

Brescia, Italy 

3 Kapitza Institute for Physical Problems, ul. Kosygina 2, 117334 Moscow 

(February 1, 2008) 

The response of a trapped Bose-Einstein condensed gas to a density perturbation generated by a 
two-photon Bragg pulse is investigated by solving the time-dependent Gross-Pitaevskii equation. We 
calculate the total momentum imparted to the condensate as a function of both the time duration 
of the pulse and the frequency difference of the two laser beams. The role of the dynamic response 
function in characterizing the time evolution of the system is pointed out, with special emphasis 
to the phonon regime. Numerical simulations are compared with the predictions of local density 
approximation. The relevance of our results for the interpretation of current experiments is also 
' discussed. 
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I. INTRODUCTION 
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Recent experiments based on two-photon Bragg spectroscopy have given first valuable information on the 
density-density response function of a trapped Bose-Einstein condensate, thereby providing an important opportunity 
to understand the role of quantum correlations and to check theoretical predictions on the dynamics of nonuniform 
condensates at both low and high excitation energy |^]. In these experiments two laser beams impinge upon a 
condensate, whose atoms can undergo stimulated light scattering events by absorbing a photon from one of the beams 
and emitting into the other. The difference in the wave vectors of the beams defines the momentum transfer in a 
single scattering event, hq , while the frequency difference defines the corresponding energy transfer, hw. Both the 
values of q and to can be tuned by changing the angle between the two beams and varying their frequency difference. 
For small values of q the system is excited in the phonon regime and the response is usefully detected by measuring 
the net momentum imparted to the gas. The latter is measured by observing the center of mass motion of the sample 
after releasing the trap. 

Since the momentum acquired by the condensate is the physical quantity measured in these experiments, it is 
important to understand what are the underlying mechanisms characterizing its behavior and what can be learnt 
from a careful study of its time dependence. In this work we investigate this problem by using the Gross-Pitaevskii 
(GP) theory for the order parameter of the condensate, with special emphasis on the phonon (low q) regime. First 
we present results of simulations based on the numerical integration of the time dependent GP equation. Then 
we analyze the same results from the viewpoint of linear response theory, testing the accuracy of the local density 
approximation for the dynamic structure factor. Finally we compare our results with the experimental data obtained 
from two-photon Bragg spectroscopy Q . 

II. TIME DEPENDENT GROSS-PITAEVSKII EQUATION 

The interaction of the atoms with the laser field is described by the Hamiltonian 

ffBra gg = ^(^e-^ t + 5p q e- i ) (1) 
where Spq = p q — (pq) is the fluctuation of the density operator 
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while the strength of the perturbation, V, is proportional to the intensity of the laser beam and can be related to the 
two- photon Rabi frequency The validity of expression (|l|) is not limited to the regime of small perturbations, but 
can also be used to explore non linear effects, including, for example, the coherent splitting of the condensate 0. 
We simulate the time evolution of the system by solving the Gross-Pitaevskii (GP) equation 



— V 2 + K (r) + g|$| 2 + 0(t)V cos{qz - ut) 
Zm 



*(r,t) 
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for the order parameter of the gas, <&(r,i). The harmonic confinement is represented by the potential Vho( r ) = 
(m/2)[u) 2 L r 2 L +ulz 2 }, with r\ = x 2 + y 2 1 while the coupling constant in the mean field term, g = Airf^a/m, is fixed by 
the s-wave scattering length a. The last term in the r.h.s. comes from Eq. (|l|), while is the Heaviside step function. 
The GP equation (0) is expected to provide an accurate description of the system at sufficiently low temperature 
We numerically solve it for a system of sodium atoms confined by a cigar-shaped potential, with co± — 8uj z . For 
the Thomas-Fermi parameter 15Na/a± we choose the value 10 4 . Here a± = y h/mu>± is the oscillator length in the 
transverse direction. Actually, by a proper rescaling of the spatial coordinates (in units of a±) and energies (in units 
of fouj±), the solutions of the time-dependent GP equation (0) depend only on the Thomas-Fermi parameter, the ratio 
<^z/<^j_ and, of course, the parameters characterizing the Bragg potential. Finally, we take q along z to preserve the 
axial symmetry of the problem as in the MIT experiment In order to select the phonon regime, one has to choose 
a wave vector q such that q£ < 1, where £ is the healing length. The latter is defined as £ = [&ima\~ 1 t 2 , where n is 
the density evaluated in the center of the trap. In our calculation, a typical value is q = SaJ 1 , with qt; = 0.42. 

The numerical integration of the time dependent GP equation (|3|) is performed by using a method previously 
developed in , which is suitable for axially symmetric condensates. The equation is mapped on a two-dimensional 
grid of points and solved by using a combination of Fast Fourier and Crank-Nicholson integration algorithms. From 
the solution of Eq. ([?]) , one can easily calculate net momentum imparted to the condensate by direct integration of 
the current density associated with the order parameter 



dr$*(r,i)V z $(r,i) + c.c. 
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In Fig. [I] we show the results obtained for P z (t), having fixed the frequency difference of the Bragg pulses to be 
lu = I3u>±. We use different values of the strength V in order to check whether we are in the appropriate limit of 
linear response to the external perturbation. In the linear regime the momentum P z of the condensate should depend 
quadratically on the intensity V of the Bragg pulse. We have verified that, by taking V smaller than Q.2Hlu±, the 
linear regime is ensured for times shorter than 0.2 (2tv/uj±). We have solved the time dependent GP equation also 
in the case of a strong Bragg pulse (V > lOhio^). For such a pulse the change in the momentum distribution of the 
condensate is huge (J] and our simulation shows that, in this case, the momentum transfer stays systematically below 
the value predicted by linear theory. 

In Fig. |2j the net momentum transfer is shown as a function of the frequency ui for two different choices of the time 
duration tb of the Bragg pulse. The figure shows that the shape of the signal exhibits a significant time dependence, 
becoming narrower by increasing tb ■ 



III. LINEAR RESPONSE THEORY 



In order to better understand the behavior of the momentum transfer as a function of t and u>, we use the formalism 
of linear response function, valid for weak intensities of the Bragg pulse. 
First, we write the Heisenberg equation of motion for the momentum P z = 

' //M/;) =i<[f>^o t ]> (5) 



dt ih 



where the Hamiltonian of the system is 
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H ^ = E h + E ^ r *) + 9 E S & ^W^Bragg • (6) 
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The commutator can be explicitly evaluated, obtaining the exact equation 



^1 = -m^Z *f {{5p\)e~^ - (Spje^) . (7) 

The first term on the r.h.s originates from the confining oscillator potential, while the second one from the Bragg 
perturbation. The quantity Z = 1S the expectation value of the 2-th component of the center of mass 

coordinate. Notice that the commutator in Eq. (^|) is not affected by the kinetic energy term, nor by two-body 
interactions which are translational invariant quantities. 

Now, one can use the linear response theory for the density fluctuations. A convenient way consists in writing the 
time dependent perturbation as 

iV f p—iw't 
tf(i)F Bragg = —8p\ / dJ — + H.c. , (8) 

4-7T M J LJ' — UJ + IT] 

which follows from the Fourier representation of the step function. This perturbation induces density fluctuations in 



the form (l0| 
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where 



<<W = -r- / du,' [5(q, u/) - S(-d, -co')} - - , (9) 
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is the dynamic structure factor of the system, with f3 — 1/kT, hu) mn — E m — E n and Z is the usual partition function. 
In deriving this result, one assumes that the fluctuation (Spq) induced by the e~ lu t component of the perturbation 
(B) oscillates like e~ llAj *, so neglecting the term oscillating like e lu The latter exactly vanishes in a uniform body 
and, in general, is exponentially small if the value of q is much larger than the inverse of the size of the system. In 
the following we will always neglect this contribution. 

Inserting the density fluctuation (||) into Eq. (Q), one gets the following equation for the rate of momentum transfer 
at t > 0: 

dPS) _ 2ry , 2g (V\ 2 f ,..wn[(u-u')t] 



This equation explicitly contains the effect of the confining potential through the first term in the r.h.s.; this term is 
small for times short compared to the oscillator period 2it/lo z . We have verified that for the times considered in our 
simulation the effect of the external potential in the equation for the momentum rate can be safely ignored. 



For very short times the expansion of (11) yields the model independent result 

dP x {t) _ LoV 2 qH 3 
dt 3m 

where we have made use of both the /-sum rule J duj lo S '(q, ui) — Nhq 2 /2m and the relationships 



(12) 



J du [S{q, lu) - S(-q, -«)] = J dtouj 2 [5(q, u) - 5(-q, -w)] = . (13) 



A closed equation for P z (t) can be obtained by taking the time derivative of (|llj) and using the exact equation 
dZ(t)/dt — P z (t)/m. If the duration of the pulse is short compared to the oscillator time, but at the same time large 
compared to the inverse of the frequency of the applied field, then the equation for the momentum rate approaches 
the golden rule result: 



dP z (t) (V\ 2 2tt 



dt \ 2 I h 



[S(q,io)-S(-q,-Lo)} . (14) 



Although restrictive, the conditions ut ^> 1 and u z t 1, needed to derive result (pT|), are compatible. In principle, 
even in the phonon regime where the excitation energy Tiw should be smaller than the chemical potential /i, the two 
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conditions can be simultaneously satisfied if one works in the Thomas- Fermi regime where /i 3> ftiu±. In practice, we 
have found difficulties in exploring the large time domain because of the occurrence of significant nonlinear effects. 

Analogous results, valid in the linear regime, can be obtained for the energy rate. In this case one finds the standard 
result of perturbation theory, 

dE{t) V 2 f sin [{cj - u')t] 

at 2n J to — to' 

which yields the golden rule result 

dE(t) (V\ 2 2iz 



dt u [l) (16) 

in the large t limit. From the experimental viewpoint it is very difficult to extract the net energy transfer since it 
would require a very high precision in the determination of the release energy. Furthermore, the release energy does 
not coincide with the total energy of the sample which includes also the contribution of the confining potential. 

It is worth noticing that result (|ll|) for the momentum rate, like Eq. ( |l5| ) for the energy rate, is sensitive to the 
difference [S(q, u>) — S(— q, — o>)] rather than to the dynamic structure factor itself. This difference characterizes the 
imaginary part of the density-density response function 

Im ( Xq H) = -(n/h) [5(q, to) - S(-q, -u)} (17) 

and follows from the fact that atoms can scatter by absorbing a photon from either of the laser beams. This represents 
an important difference with respect to other scattering experiments (like, for example, neutron scattering from helium) 
where, by detecting the scattered probe, one instead measures directly the dynamic structure factor. The dynamic 
structure factor and the imaginary part of X q (w) have the same behavior at T = since in this case 5*(q, lj) is zero 
for negative to. They instead differ at finite temperatures if fc^T is of the order or higher than the excitation energy 
hu>. Actually the difference ( |l7| ) significantly suppresses the thermal effects exhibited by the dynamic structure factor 
so that, by measuring x q (w), one has an easy access to the zero temperature value of 5*(q, to). For this reason, even if 
experiments are carried out at temperatures which do not satisfy the condition kT <C fuv, one can safely restrict the 
theoretical analysis to the simpler T — case, provided the temperature is sufficiently low to ignore the effects due 
to the thermal depletion of the condensate. 



IV. LOCAL DENSITY APPROXIMATION 



If the value of the wave vector q is larger than the inverse of the size of the condensate the response of the system 
can be obtained, in first approximation, using a local density approximation (LDA). In fact, in this case the effects 
of discretization in the excitation spectrum can be safely ignored and the system can be treated as a locally uniform 
medium. In our simulation we have qR± — 21 and qR z — 168 and consequently this approximation should be well 
satisfied. In the local density approximation the dynamic structure factor is given by the analytic expression 



<SLDA(q, w ) = 



15fi (h 2 m 2 - E 2 T ) 
"13 E r [i 2 



1 - 



(H 2 lu 2 - E 2 ) 



-,1/2 



2E r [i 



(18) 



where 
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E r 



K 2 q 
2m 



(19) 



is the free recoil energy and /i = ft 2 /(2m!; 2 ) is the Thomas-Fermi value of the chemical potential. In deriving result 
( |l8| ) one has to average the Bogoliubov expression S(q,ui) — [h 2 q 2 /(2me(q))]5(ui — e(q)/h) over the Thomas-Fermi 
density profile n(r) — (l/g)[/i — T4 xt (r)], by evaluating the Bogoliubov excitation spectrum 



2m V 2m 



K l q 2 fh 2 q 2 



gn(r] 



1/2 



(20) 



at the corresponding density. In order to apply the LDA, the momentum transfer hq should not be however too large, 
because this approximation ignores the Doppler effect associated with the spreading of the momentum distribution of 
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the condensate which is expected to become the leading effect at very large values of q ||. This happens for values 
of q much larger than the inverse of the healing length, a situation that is not considered in the present work. 

Equation (18) shows that, differently from the case of a uniform gas, the dynamic structure factor of a trapped 
condensate is no longer a delta function, its value being different from zero in the interval E r < E < E r [l + 2/i/ 'Er] 1 / 2 . 
The value E = E r corresponds to the excitation energy in the region near the border where the gas is extremely 
dilute and hence noninteracting. The value E — E r [l + 2fi/ 'Er] 1 ' 2 is instead the excitation energy of a Bogoliubov gas 
evaluated at the central density. Notice that the LDA expression ( |l8| ) for S(q, to) does not depend on the direction of 
the vector q even in the presence of a deformed trap. 

Results (|llj) and ( |l8| ) allows one to evaluate the time dependence of the momentum imparted to the condensate in 
the linear regime. In Figs. and | we compare the results of this local density approximation with the ones of the 
numerical integration of the GP equation. The agreement is excellent and proves that the LDA is quite adequate to 
describe the response of the system in the conditions considered in our analysis. We have also carried out the same 
comparison for a smaller value of momentum transfer (q = la^ 1 ) finding a similar good agreement. 

The successful comparison between the predictions of the LDA and the numerical simulation suggests that equations 
( |ilt ) and ( fill) are indeed useful for the analysis of the experimental results of Ref ||] (see Fig. The values of the 
parameters used in the experiment (q = 4.3a]! 1 , u> z /lj± = 0.12, q£ — 0.39) do not differ significantly from the ones used 
in our simulation, except for the value of the chemical potential which is significantly larger (/i = 60huj±). Because 
of this, the numerical simulation of the MIT experiment would be a highly time consuming calculation, beyond the 
purpose of this work. On the other hand, in the above experimental conditions the applicability of the LDA should 
be even better, because of the increase of the value of qR± and of qR z . 

In Fig. H we compare the LDA prediction (solid line) with the available experimental data (points) . The latter are 
given in arbitrary units, since the value of V used in Q is not available, and hence the comparison is limited to the 
shape and position of the peak. The LDA curve is obtained by integrating Eq. (|ll|) between t = and t — tb and 
using expression ( |l8| ) for the dynamic structure factor. The agreement is reasonably good. It is also interesting to 
compare the LDA prediction with the golden rule (|lj) integrated over the same time interval (dashed line). The LDA 
curve is broadened and quenched as a result of the finite duration of the Bragg pulse, which enters Eq. (y_lj) through 
the sinusoidal factor in the integrand. In the golden rule (Q) the same factor is replaced with a delta function, but 
the figure shows that the effect of this replacement is significant and can not be neglected. 



V. CONCLUSIONS 



In conclusion we have shown that the net momentum imparted to a trapped condensate by light scattering can 
be calculated by solving the time dependent Gross-Pitaevskii equation in a regime which is significant for current 
experiments. In the linear response limit, this quantity is directly connected with the dynamic structure factor for 
which the local density approximation (LDA) turns out to be a reliable approximation in the explored range of <?'s. 
The connection between the dynamic structure factor and the quantities measured in stimulated light scattering 
experiments in the diverse accessible regimes deserves further investigations (see also Ref. ||ll|). The present work 
shows that, within appropriate limits, a quantitative analysis of experimental data is indeed possible. More accurate 
and systematic measurements, including the explicit dependence of the momentum imparted to the condensate on 
the duration of the pulse, would allow for a better test of the theoretical predictions for the dynamic structure factor 
of trapped Bose-Einstein condensates. 
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FIG. 1. The quantity (A/V 2 )P Z , in units of [ftw^ax] -1 i s plotted against time, in units of vj 1 — 2ty/lu±. The Bragg pulse 
is applied starting at t = 0, with q — 3a]] 1 and uj = 13o;x. The dot-dashed line is the model independent short time expansion 
(fl2|). The solid line is the prediction of LDA, obtained by using the dynamic structure factor ( |l8| ) in Eq. ([H). The two dashed 
lines are the results of the numerical integration of the GP equation ([]), with two different values of the Bragg intensity V , 
namely, V = 0.2hco± (long dashed) and 0Ahuj± (short-dashed). 
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FIG. 2. The quantity (4/V 2 )P z , in units of [hu} 2 ± a ± ]' 1 is plotted against the frequency difference of the two laser beams, uj 
in units of ujj_, for q — SaJ 1 and two different values of the time duration of the Bragg pulse tb, in units of vj_ = 2tt/ui±. Solid 
lines are the predictions of LDA, as in Fig. 1. Dashed lines are the results of the numerical integration of the GP equation (|^), 
with Bragg intensity V = 0.2htu±. 




u/2n [kHz] 

FIG. 3. The quantity (4/V 2 )P 2 , in units of [hu 2 ± ax] 1 is plotted against frequency, lu/2ty, for Q__= 4.3a x 1 and 
tb ~ 0.06(27r/cjj_) = 400 ^tsec. The solid line is the prediction of LDA, obtained from Eqs. (|llj) and (O) using the pa- 
rameters of Ref. jj. The dotted line is the prediction of the golden rule (Q). The experimental values of Ref^| (filled circles) 
are given in arbitrary units. 
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